#-*- coding: utf-8 -*- 

'''
Created on Feb 16, 2011
k Means Clustering for Ch10 of Machine Learning in Action
@author: Peter Harrington
'''
from numpy import *
import matplotlib.pyplot as plt

#读数据
def loadDataSet(fileName):        
    dataMat = []   #创建列表，存储读取的数据
    fr = open(fileName)
    for line in fr.readlines(): #读每一行
        line1=line.strip()     #删头尾空白
        curLine = line1.split('\t') #以\t为分割，返回一个list列表
        fltLine = map(float,curLine)#str 转成  float
        dataMat.append(fltLine)     #将元素添加到列表尾
    return dataMat 

#算距离
def distEclud(vecA, vecB):                  #两个向量间欧式距离
    return sqrt(sum(power(vecA - vecB, 2))) #la.norm(vecA-vecB)


#初始化聚类中心：随机产生区间范围的值作为新的中心点
def randCent(dataSet, k):
    #特征维度
    n = shape(dataSet)[1]   
    #创建聚类中心的矩阵 k x n                
    centroids = mat(zeros((k, n)))
    #遍历n维特征         
    for j in range(n):     
        #第j维特征属性值min   ,1x1矩阵                 
        minJ = min(dataSet[:, j])
        #区间值max-min，float数值    
        rangeJ = float(max(dataSet[:, j]) - minJ)
        #第j维，每次随机生成k个中心 
        centroids[:, j] = mat(minJ + rangeJ * random.rand(k, 1))
    return centroids

#随机抽取两个样本，初始的中心质点
def randChosenCent(dataSet,k):
    m=shape(dataSet)[0] #样本数
    centroidsIndex=[]
    dataIndex=range(m)#样本的序号
    for i in range(k):
        randIndex=int(random.uniform(0,len(dataIndex)))
        centroidsIndex.append(dataIndex[randIndex])
        #删除已经被抽中的样本
        del(dataIndex[randIndex])
    centroids=dataSet[centroidsIndex]
    return mat(centroids)


	
#kmeans算法:当所有样本的分配结果不发生改变时，终止迭代
def kMeans(dataSet, k, distMeas=distEclud, createCent=randChosenCent):
    m = shape(dataSet)[0]   #样本总数
    #分配样本到最近的簇：存[簇序号,距离的平方] 
    clusterAssment = mat(zeros((m,2)))       
    #step1:#初始化聚类中心                                     
    centroids = createCent(dataSet, k)   
    clusterChanged = True
    #所有样本分配结果不再改变，迭代终止
    while clusterChanged:   
        clusterChanged = False        
        #step2:分配到最近的聚类中心对应的簇中
        for i in range(m):   
            minDist = inf; minIndex = -1  #对于每个样本，定义最小距离
            for j in range(k):  #计算每个样本与k个中心点距离
                distJI = distMeas(centroids[j,:],dataSet[i,:]) 
                if distJI < minDist: 
                    minDist = distJI; minIndex = j  #获取最小距离,及对应的簇序号
            if clusterAssment[i,0] != minIndex: clusterChanged = True 
            clusterAssment[i,:] = minIndex,minDist**2 #分配样本到最近的簇
        print 'centroids=',centroids
        # sse=sum(clusterAssment[:,1])
        # print 'sse=',sse
        #step3:更新聚类中心
        for cent in range(k):#样本分配结束后，重新计算聚类中心
            #获取该簇所有的样本点
            ptsInClust = dataSet[nonzero(clusterAssment[:,0].A==cent)[0]]
            #更新聚类中心：mean()对array类型求均值，axis=0沿列方向求均值。
            centroids[cent,:] = mean(ptsInClust, axis=0) 
    return centroids, clusterAssment


#当误差平方和SSE收敛，退出循环
def kMeansSSE(dataSet, k, distMeas=distEclud, createCent=randChosenCent):
    m = shape(dataSet)[0]
    #分配样本到最近的簇：存[簇序号, 距离的平方]
    clusterAssment=mat(zeros((m,2)))
    #step1:#初始化聚类中心
    centroids = createCent(dataSet, k)
    sseOld=0
    sseNew=inf
    while(abs(sseNew-sseOld)>0.001):
        sseOld=sseNew
        #step2:将样本分配到最近的质心对应的簇中
        for i in range(m):
            minDist=inf; minIndex=-1
            for j in range(k):
                #计算第i个样本与第j个质心之间的距离
                distJI = distMeas(centroids[j,:], dataSet[i,:])
                #获取到第i样本最近的质心的距离,及对应簇序号
                if distJI<minDist:
                    minDist=distJI; minIndex=j
            clusterAssment[i,:] = minIndex, minDist**2 #分配样本到最近的簇
        print 'centroids=',centroids
        sseNew = sum(clusterAssment[:, 1])
        print 'sse=', sseNew
        #step3:更新聚类中心
        for cent in range(k):#样本分配结束后，重新计算聚类中心
            ptsInClust=dataSet[nonzero(clusterAssment[:,0].A==cent)[0]]
            #按列取平均,mean()对array类型
            centroids[cent,:] = mean(ptsInClust, axis=0)
    return centroids, clusterAssment



#2维数据聚类效果显示
def datashow(dataSet,k,centroids,clusterAssment):  #二维空间显示聚类结果
    from matplotlib import pyplot as plt
    num,dim=shape(dataSet)  #样本数num ,维数dim
    
    if dim!=2:
        print 'sorry,the dimension of your dataset is not 2!'
        return 1
    marksamples=['or','ob','og','ok','^r','sb','<g'] #样本图形标记
    if k>len(marksamples):
        print 'sorry,your k is too large,please add length of the marksample!'
        return 1
    
    #绘所有样本
    for i in range(num):
        markindex=int(clusterAssment[i,0])#矩阵形式转为int值, 簇序号
        #特征维对应坐标轴x,y；样本图形标记及大小
        plt.plot(dataSet[i,0],dataSet[i,1],marksamples[markindex],markersize=6)
    #绘中心点            
    markcentroids=['dr','db','dg','dk','^b','sk','<r']#聚类中心图形标记
    for i in range(k):
        plt.plot(centroids[i,0],centroids[i,1],markcentroids[i],markersize=15)
   
    plt.title('k-means cluster result') #标题        
    plt.show()

#聚类前，绘制原始的样本点
def originalDatashow(dataSet):
        #样本的个数和特征维数
    num,dim=shape(dataSet)
    marksamples=['ob'] #样本图形标记
    for i in range(num):
        plt.plot(datamat[i,0],datamat[i,1],marksamples[0],markersize=5)
    plt.title('original dataset') #标题
    plt.show()

if __name__=='__main__':
    #=====kmeans聚类
    # #获取样本数据
    datamat=mat(loadDataSet('testSet.txt'))
    #原始数据显示
    originalDatashow(datamat)

    #*****kmeans聚类
    k=4 #用户定义聚类数
    mycentroids,clusterAssment = kMeansSSE(datamat, k)
    # #绘图显示
    datashow(datamat, k, mycentroids, clusterAssment)
		


